/*
 * util.c
 *
 *  Created on: Oct 25, 2013
 *      Author: jiang
 */

#include <stdio.h>
#include "bntseq.h"
#include "kseq.h"
#include "util.h"
#include "global.h"

//DEBUG
#include <inttypes.h>
//DEBUG

#define _set_pac(pac, l, c) ((pac)[(l)>>2] |= (c)<<((~(l)&3)<<1))
#define _get_pac(pac, l) ((pac)[(l)>>2]>>((~(l)&3)<<1)&3)

uint8_t *add1(const kseq_t *seq, bntseq_t *bns, uint8_t *pac, int64_t *m_pac, int *m_seqs, int *m_holes, bntamb1_t **q)
{
	bntann1_t *p;
	int i, lasts;
	if (bns->n_seqs == *m_seqs) {
		*m_seqs <<= 1;
		bns->anns = (bntann1_t*)realloc(bns->anns, *m_seqs * sizeof(bntann1_t));
	}
	p = bns->anns + bns->n_seqs;
	p->name = strdup((char*)seq->name.s);
	p->anno = seq->comment.s? strdup((char*)seq->comment.s) : strdup("(null)");
	p->gi = 0; p->len = seq->seq.l;
	p->offset = (bns->n_seqs == 0)? 0 : (p-1)->offset + (p-1)->len;
	p->n_ambs = 0;
	for (i = lasts = 0; i < seq->seq.l; ++i) {
		int c = nst_nt4_table[(int)seq->seq.s[i]];
		if (c >= 4) { // N
			if (lasts == seq->seq.s[i]) { // contiguous N
				++(*q)->len;
			} else {
				if (bns->n_holes == *m_holes) {
					(*m_holes) <<= 1;
					bns->ambs = (bntamb1_t*)realloc(bns->ambs, (*m_holes) * sizeof(bntamb1_t));
				}
				*q = bns->ambs + bns->n_holes;
				(*q)->len = 1;
				(*q)->offset = p->offset + i;
				(*q)->amb = seq->seq.s[i];
				++p->n_ambs;
				++bns->n_holes;
			}
		}
		lasts = seq->seq.s[i];
		{ // fill buffer
			if (c >= 4) c = lrand48()&3;
			if (bns->l_pac == *m_pac) { // double the pac size
				*m_pac <<= 1;
				pac = realloc(pac, *m_pac/4);
				memset(pac + bns->l_pac/4, 0, (*m_pac - bns->l_pac)/4);
			}
			_set_pac(pac, bns->l_pac, c);
			++bns->l_pac;
		}
	}
	++bns->n_seqs;
	return pac;
}


int getKmerStrandIndicatorBit(uint32_t kmerMapArrayPos,const uint8_t *kmerStrandIndicator)
{
	uint32_t arrayIndex = kmerMapArrayPos/8;
	int shift = kmerMapArrayPos%8; //ZERO based, from left to right
	uint8_t number = kmerStrandIndicator[arrayIndex];
	int bit = number & (1 << (7 - shift));
	return bit;
}
int setKmerStrandIndicatorBit(uint32_t kmerMapArrayPos,uint8_t *kmerStrandIndicator)
{
	uint32_t arrayIndex = kmerMapArrayPos/8;
	int shift = kmerMapArrayPos%8; //ZERO based, from left to right

	//fprintf(stderr,"kmerMapArrayPos: %llu, shift: %d ", kmerMapArrayPos,shift);
	//fprintf(stderr,"Before: %x ", kmerStrandIndicator[arrayIndex]);
	kmerStrandIndicator[arrayIndex]  |= 1 << (7 - shift);
	//fprintf(stderr,"After: %x ", kmerStrandIndicator[arrayIndex]);
	//fprintf(stderr,"check bit: %d \n", getKmerStrandIndicatorBit(kmerMapArrayPos,kmerStrandIndicator));
	//fprintf(stderr,"check bit one more: %d \n", getKmerStrandIndicatorBit(kmerMapArrayPos-1,kmerStrandIndicator));
	return 0;
}


uint8_t *getTarget(const uint8_t *pac, int64_t l_ref_oneStrand, uint64_t position, uint32_t targetLen)
{
	uint64_t i,pos, leftPos_PlusStrand;
	uint8_t *buf, *rvc_buf;

	buf = (uint8_t*)calloc(targetLen + 1, 1);
	rvc_buf = (uint8_t*)calloc(targetLen + 1, 1);


	//fprintf(stderr,"%" PRId64 "\n",position);
	//fprintf(stderr,"%" PRId64 "\n",position -BIGNUMBER);
	//

	if (position < BIGNUMBER)
	{

		if (position < l_ref_oneStrand *2 )
		{
			//fprintf(stderr,"+ working \n");
			//fprintf(stderr,"pos: %" PRId64 "\ntarget:",position);
			for (i = 0; i < targetLen; ++i) {
						pos = i + position;
						buf[i] = pac[pos>>2] >> ((3 - (pos&3)) << 1) & 3;
						//fprintf(stderr,"%d",buf[i]);
			}
		}

		else{
			for (i = 0; i < targetLen; ++i) {
				buf[i] = 1;
				//fprintf(stderr,"%d",buf[i]);
			}
		}
			free(rvc_buf);
			return buf;


	}
	else{ //the map pos located in - strand. Do reverse complement for the target


		leftPos_PlusStrand = (l_ref_oneStrand -1) - (position-BIGNUMBER) - (targetLen -1);
		if (leftPos_PlusStrand < l_ref_oneStrand *2 ){
		//DEBUG
		//fprintf(stdout,"rev_pos: %llu",leftPos_PlusStrand);
		for (i = 0; i < targetLen; ++i) {
				pos = i + leftPos_PlusStrand;
				buf[i] = pac[pos>>2] >> ((3 - (pos&3)) << 1) & 3;
		}
		// reverse complement the buf
		for (i = 0; i < targetLen; ++i)
				rvc_buf[targetLen - i - 1] = 3 - buf[i];
		}
		else{
			for (i = 0; i < targetLen; ++i) {
				rvc_buf[i] = 1;
							//fprintf(stderr,"%d",buf[i]);
			}
		}
		free(buf);
		return rvc_buf;


	}

}

uint32_t hammingDistance(uint8_t *a, uint8_t *b,  uint32_t len)
{
	uint32_t i, dist =0;
	//fprintf(stderr,"\nCOMPARE: ");
	for( i=0; i< len; i++)
	{
		if((uint8_t) a[i] != (uint8_t)b[i])
		{   //fprintf(stderr," %d!= %d ",a[i],b[i]);
			dist++;

		}

	}
	return dist;
}



int compareParSecondDec(const pair *s1, const pair *s2)
{

	 if (s1->second < s2->second) return 1;
	 if (s1->second ==  s2->second) return 0;
	 if (s1->second > s2->second) return -1;

}
